Limits on deviations from the inverse-square law on megaparsec scales 
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We present an attempt to constrain deviations from the gravitational inverse-square law on 
large-scale structure scales. A perturbed law modifies the Poisson equation, which implies a scale- 
dependent growth of overdensities in the linear regime and thus modifies the power spectrum shape. 
We use two large-scale structure surveys (the Sloan Digital Sky survey and the Anglo- Australian 
Two-degree field galaxy redshift survey) to constrain the parameters of two simple modifications of 
the inverse-square law. We find no evidence for deviations from normal gravity on the scales probed 
by these surveys (~ 10 23 m.) 
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INTRODUCTION 

Einstein's theory of general relativity, which general- 
izes Newton's law of gravity, has been extensively tested 
and verified. But precision tests of the inverse-square 
force law have been performed only on scales < 10 13 
m (P| and references therein). An enormous extrapo- 
lation is required to apply this law to cosmological scales 
> 10 22 m 2]. 

Some recent theories propose a deviation from New- 
tonian gravity at large distances, perhaps due to ex- 
tra dimensions, branes, or new particles (e.g.jjj)- These 
theories, mostly motivated by finding an alternative ex- 
planation to dark energy for the cosmic acceleration 
(e.g.0, IE IE 0)> modify gravity on scales comparable 
to the horizon scale ~ 10 26 m. 

Here we set out to constrain theories in which gravity 
deviates from the inverse-square law on scales of ~10 Mpc 
or 10 23 m. We evolve the power spectrum (and bispec- 
trum) under a perturbed law of gravity in the linear (and 
mildly non-linear) regimes, from recombination to the 
present day. The perturbed law modifies the evolution 
of density fluctuations, altering, in particular, the power 
spectrum shape. We compare the prediction for the mod- 
ified power spectrum to that measured from two galaxy 
redshift surveys: the Anglo-Australian Two-degree Field 
Galaxy Redshift Survey (2dFGRS;Jl@) and the Sloan 
Digital Sky Survey (SDSS; [IE El El ) . 

In the models we consider, gravity is consistent with 
the inverse-square law on scales smaller than ~1 Mpc, 
so the deviation does not affect early universe physics 
through recombination. Thus, we assume that the stan- 
dard computation of the cosmic microwave background 
(CMB) anisotropies holds, and that the CMB obser- 
vations provide the mass power spectrum (and bispec- 
trum) at recombination. Moreover, as explained later, 
we marginalize over reasonable priors for the power 
spectrum amplitude and the primordial power spectrum 
slope, thus making our analysis insensitive to reasonable 
changes in the other cosmological parameters. We expect 
deviations from the inverse-square law to modify the inte- 



grated Sachs- Wolfe (ISW) effect. However, as this affects 
the CMB on large scales where cosmic variance is large, 
we do not discuss ISW here. 

In the context of Newtonian gravity, a similar ap- 
proach has been widely used in the literature to con- 
strain cosmological parameters (e.g., EH) and biasing, 
the relation between clustering of galaxies and cluster- 
ing of mass (e.g.. [13. Hfij|). Here, we fix the background 
cosmology E3 and assume a scale- independent bias [IE] • 
The deviation from Newtonian gravity is parameterized 
by a length-scale and a strength on which we place con- 
straints. We find no compelling evidence for deviations 
from the inverse-square force law on large-scale structure 
scales. This null result leaves open the possibility that 
both Newton's gravity and our weak assumptions about 
scale-independent bias and the primordial power spec- 
trum are incorrect in some contrived way so that their 
effects conspire to exactly cancel. We regard such a sce- 
nario as exceedingly unlikely. This is the first attempt to 
use cosmological data to constrain gravity on these scales, 
but for related theoretical work, see [ill HH Hsl lil| . 



METHOD 

We consider a force of gravity that deviates from an 
inverse-square law by a small perturbation on Mpc scales. 
We examine two possible functional forms for this per- 
turbation, each parameterized by a strength and a length 
scale. Since Newtonian gravity has been tested on very 
small (mm) scales up to solar system scales Q , both mod- 
els recover the inverse-square law at small scales. The 
expansion rate a(t) is set by WMAP's best-fit ACDM 
background cosmology EH- Note that these models dif- 
fer from the modification to general relativity considered 
in [2fJ. Since we are testing a null hypothesis, the spe- 
cific form of the parameterization does not matter; our 
choices are based loosely on parameterizations used in 
small-scale tests of gravity. 

In the first model, we consider a Yukawa-like contribu- 
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tion to the potential, 



$(f) = -G J d 3 r' 
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where a is small. The Yukawa parameterization has been 
used in small-scale tests of gravity, where it is physically 
motivated by the exchange of virtual bosons [lj. This 
model nicely behaves like an inverse-square law at very 
large and very small length scales relative to the length 
scale A. 

In the second model (hereafter the "PL model"), we 
consider a power-law-like (PL) potential of the form, 



$(f) = -G J d 3 r 



, P(r') 



1 + e(|r — 



R 



N-l 



(2) 

where efr) = 6(r — R) is a step function, e <C 1, R is a 
physical length scale, e is the strength of the deviation, 
and we discuss N = 2. The r dependence of e ensures 
agreement with standard gravity at small scales. Inciden- 
tally, this form of modification (with constant e) can be 
generated by the simultaneous exchange of two massless 
scalar particles (pj and references therein). However, 
such effects would likely occur at much smaller scales 
than we are probing. 

Let's consider the first model. In comoving coordi- 
nates, x = ^y, the comoving gravitational potential is 
given by (f>(x) — — Gpa A.(x), where 
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, p is the average density of the universe, 



and we scale a(t) so a = 1 today. 

Transforming A into Fourier space yields 
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with F Al (s) = j^t. 

Following the same analysis as the first model, the PL 
model yields 



A(fc) = 



47T(5(fc) 

k 2 



1 + (F ^(m) 



(5) 



with Fa 2 (s) = - \% — Si , where Si(ir) denotes the 
sine integral. 



Linear theory 

In this section we compute the density fluctuations 5, 
the comoving velocity field, and the power spectrum for 
the Yukawa-like model. The results for the PL model are 



obtained simply by replacing a with e, A with R, and 
F Al (s) with Fa 2 (s). 

The universe's expansion rate is dominated by physics 
at horizon scales, beyond the scales we consider here. We 
thus use Friedman equations for the background cosmol- 
ogy. Although for the PL model we simply make this 
assumption, for the Yukawa-like model this assumption 
can be justified. The Yukawa-like model recovers the 
inverse-square law on horizon scales. A change in the 
value of Newton's constant G, as in the large-scale limit 
of the Yukawa-like model, has the effect of changing the 
overall amplitude of the growth factor, and thus of the 
power spectrum normalization, over which we marginal- 
ize. We are sensitive only to changes in the shape of 
the power spectrum on supercluser scales, which is most 
strongly affected by the growth of perturbations against 
the background cosmology, and not by the background 
cosmology itself. 

The equation for 8 to first order in perturbation theory 

is 

5 + 2-5-Gpk 2 A = 0. (6) 
a 

For a = the equation is separable, with solution}^ 

S A (k,t) = A(t)S (k). (7) 
We now look for a solution of the form 

6(k, t) = 6 A (k, t) \l + ad(k, t) (8) 

to first order in a, and solve for d(k, a(t)) The initial 
conditions d(Q) — and <i'(0) = ensure that gravity 
behaves normally at small scales and that 5 is finite at 
a = 0. For the PL model, all solutions are undefined at 
a = 0. We set d and d' to zero at a sufficiently small 
scale that changing that scale by an order of magnitude 
changes d at a = 1 by significantly less than our error 
bars. 

Figure^) shows rf(tv) for the Yukawa- like model 
(dashed line) and d(j^) for the PL model (solid line) 
for a ACDM universe at a = 1. Since a positive value 
of a makes gravity stronger, we expect more clustering 
(d > 0). The step function in the PL model's poten- 
tial causes some ringing in the Fourier transform, so the 
change in the comoving potential and d both oscillate as 
functions of k. 

We find that J7a weakens the effect of a or e at large 
scales. 

We can now compute the comoving velocity to first 
order in perturbation theory using the continuity equa- 
tion and the equation of motion in Fourier space [2l|. 
For simplicity, we'll work in an Einstein-de Sitter (EdS) 
universe. We only use the velocity to compute the bis- 
pectrum, and the bispectrum is robust to changes in cos- 
mology (see, e.g., Figure (I) in [2^|1. To first order in a, 
we find 
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FIG. 1: d(-jry) and d(-^) for the Yukawa-like and PL models, 
respectively, at a = 1 for a ACDM cosmological model (f2 m = 
0.27 and ft A = 0.73) 



small values of s -and the same \a\ value- the error is less 
than a percent. We are only concerned with effects that 
change the result by orders of magnitude, not factors of 
a few, since gravity has never been tested on these scales 
before. Thus, we believe the Taylor expansion in a and 
e is sufficient for our first search for constraints. 



Mildly non-linear theory 

The equation for 5 to second order in perturbation the- 
ory is 



5 + 2-5 + ^ci> = —[hi *k i 6 + 8* k 2 4> + kikAtf * &)], 
a a 2 a z 

(15) 

where we use the Einstein summation convention over 
spacial indices The symbol "*" denotes a convolution 
in /c-space, 
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where 

F v (s) = d(s) + sd'(s). 
The power spectrum is given by 

(2Tr) 3 P(k,t)6 D (k + k r ) = (6{k,t)S(k',t)), 

where 

(\6(k,t)\ 2 ) = (\6 A (k,t)(l + ad(k,t))\ 2 ) 



(9) 



(10) 



(11) 



(12) 



and () denote an ensemble average over all realizations 
of the universe. The present-day power spectrum derived 
using our model of modified gravity, to first order in a, 
is then 



P(k) = P Q {k) 



l + 2ad 



1 

kX 



(13) 



where P is the power spectrum expected for Newtonian 
gravity. 

To test the validity of our first-order approximation in 
a, we also compared d for an EdS universe to the Yukawa- 
like model's exact solution. The growing EdS solution for 
the density fluctuation is a hypergeometric function, 



6(s) = 2 F 1 



5 - V25+24a 5 + \/25 + 24a 9 
8 ' 8 'I 5 



;-s 2 ) sS (k), 
(14) 



where s = -j^. We find that the higher-order a correc 
tions contribute a small fraction most values of s 



k\ 



in our analysis. For example, for s = 100 (the largest 
s in our analysis), we get ~ 20% error for \a\ ~ 0.2; for 



f{k)*g{k) = J d'H j d i m8 u {l + m-k)f{l)g{m), (16) 

where S D (k) is the three-dimensional Dirac delta func- 
tion. We use the solution for v from first-order pertur- 
bation theory. 

Let 5^ represent the solution to 5 in linear theory, 
and 5^ represent the second-order term (of order <5|), 
so 5 = 5^ + plus smaller, higher-order terms. 

For a = (or e = 0), the solution for the second-order 
term is 



5f = A\t)8f\k) 



(17) 



= A \^-5 * 5 a + hd * -S + j-j^-S * -^5 j . 

We then follow the same approach as in the linear the- 
ory. Using the first model as our example, we look for a 
solution to first order in a of the form 



P>=*i 2) [l + a$(M)], 



(18) 



Switching the time variable to a, the resulting differential 
equation for g(k, a) in an EdS universe is 

d 2 g 11 dg 7 „ 13 A - , 

^ + Ta£ + ^=W^ N(k ' a) ' (19) 

where 

N(k,a) = M^)^ 
hi 
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FIG. 2: Results from the first (Yukawa-like) model. Left: Likelihood as a function of a and A using data up to k ~ 0.15 
h/Mpc using 2dF data. Contours denote one sigma marginalized (solid line), one sigma joint (dashed line), and two sigma 
joint (dashed-dotted line). Middle: Chi square as function of a, marginalized over A, using data up to k ~ 0.15. Right: Values 
of a corresponding to the maximum chi square with error bars denoting one sigma, as a function of the maximum k included 
in the analysis. 
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with initial conditions g(k, 0) = and J^g(k, 0) = 0. The 
bispectrum is given by 

B{h, ha) 6 D (ki + k 2 + k 3 ) = {5(ki)8(k2)S{k 3 )) 

= (fc x ) + 5( 2 ) (fc 2 )] (k 2 ) + 8<n (k 2 )} 

x[8^(k 3 ) + 8^(k 3 )]}. (21) 

The leading term, (S^- 1 ' (ki)8^ (k 2 )S^ (k 3 )) , vanishes un- 
der Gaussian initial conditions. Modified gravity still 
preserves gaussianity in the linear regime, and introduces 
a correction to the bispectrum in second-order perturba- 
tion theory. To first non-vanishing order in 8 a and a, we 
get 



(S(ki)8(k 2 )S{k 3 )) = A^SoikMh^ih)) 

+ cyc. (22) 



where again we can pull d and g out of the ensemble aver- 
age, k 3 — —ki — k 2 , and eye. denotes cyclic permutations 
of the subscripts 1,2,3. 

Since we already know d, we need to solve only 
for the last term. To do so, let's examine the term 
X = (8 (k 1 )S o (k 2 )Si 2 \k 3 )g(k 3> a)). FromEq.(HHJ) we can 
write, 

d 2 X UdX 7 v /f/ r^f,r^ 3 Ar , r 
ft? + i£ Tfc + 5?* = ^(kMk 2 )^N(k 3 ,a)). 



N(k,a) is a sum of terms with k-space convolutions and 
a dependence. We can thus rewrite Eq. (|20|l as 

N{k, a) = J2 So{k)Jn(k, a) * 8 (k)K n (k, a), (24) 



where the expressions for J n and K n are easily obtained 
by equating Ea.pijl with Ea.l(20|). Using Eq.lO, we 
rewrite the right hand side of Eq. (|23|l as 



- <J d lHmH° { l + rn - k 3 )Uk)UW)8<>{™) 

(25) 



which becomes 



x J2 n Jn{l,a)K n (m,a)) 



n 

(26) 

We can then solve Ea. H23(l for each function of a that 
appears in the above sum, and get an expression for the 
bispectrum. 

The full expression for the bispectrum given by, 

S(fci,fc 2 )<5 D (fci + fc 2 + fc 3 ) - J{k u k 2 )Po{ki)Po{k 2 ) + cyc. , 

(27) 



where 



J(ki,k 2 ) = J (fci, k 2 ) + a (j (fci, k 2 ) (^(j^jj + 




(23) 
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FIG. 3: Results from the PL model. Left: Likelihood as a function of e and R using data up to k ~ 0.15 h/Mpc using 
SDSS data. Inner contours denote one sigma marginalized (~ A% 2 = 1 solid line), bottom contour denotes one sigma joint 
(~ Ax 2 = 2.3 dashed line). Right: Chi square as function of e, marginalized over R, using data up to k ~ 0.15. 




Since the bispectrum contains complementary infor- 
mation from the power spectrum, we anticipate that 
combining the two should produce stronger constraints. 
However, we leave the bispectrum analysis to future 
work. 



(28) 



J n (fci , fc?)P n (fci )P n (k?) is the Newtonian gravity result 
given by Eq.(40) of [23], and 
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Here, 2 ^ i ' u ' u,u '' denotes the partial derivative of the hy- 
pergeometric function with respect to the first argument, 
and 2 ^7' i ( ' 1 ' ' ) denotes the partial derivative with respect 
to the second argument. 



ANALYSIS 

We can use the power spectrum measurements from 
the 2dFGRS 0] and the SDSS [13] to constrain the pa- 
rameters of our models. The bispectrum derivation shows 
that the bispectrum should help disentangle the effects 
of biasing from possible modifications to gravity, but the 
analysis will follow in a subsequent paper. 

The z — power spectrum for normal gravity can be 
theoretically calculated using the transfer functions com- 
puted with CMBFAST for the WMAP best-fit parame- 
ters. The theoretical power spectrum for given values of 
a and A (or e and R) can then be obtained from Eq. (1 1 31) . 
Assuming Gaussian likelihood, we computed the likeli- 
hood for each survey's power spectrum measurements as 
a function of a and A (or e and R) . 

To make our analysis insensitive to assumptions about 
the background cosmology, we marginalized over reason- 
able priors for the power spectrum amplitude and the 
primordial power-law power spectrum slope, n. For the 
amplitude, we find that the ranges 0.6 to 1.4 (for the 
first model) and 0.4 to 1.4 (for the PL model) times the 
best-fit power spectrum amplitude for normal gravity are 
sufficient to allow for changes in the expansion rate of 
the universe. The values of f2 TO , and the dark energy 
equation of state parameter (assumed to be constant) af- 
fect primarily the power spectrum amplitude in linear 
theory and for the scales of interest. We marginalized 
over n using the WMAP prior (n = 0.99 ± 0.04 at the 
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one-sigma level |l3(). The large scale structure power 
spectrum shape has a weak dependence on T ~ Q m h, 
which shows up only on the largest scales and is strongly 
degenerate with O haryon and n 0, 0] . The large-scale 
structure power spectrum has a very weak sensitivity to 
^baryon via two effects. One is the appearance of the so- 
called baryonic wiggles; as this effect has not yet been 
detected from these data sets, it cannot yet give sensi- 
tivity to Qbaryon- The other effect is a weak change in 
power spectrum slope, degenerate with n. We are also 
somewhat insensitive to the uncertainty in neutrino mass. 
Increasing neutrino mass affects first the amplitude and 
then the shape of the galaxy power spectrum in the scales 
of interest (see, e.g., figure(8) of |25|). Present constraints 
would be weakened if a running of P(k) and/or a redshift 
dependence of the dark-energy equation of state were ex- 
plored. 

Lastly, since we are mainly interested in a constraint 
on a and/or e rather than on the scale A and/or R, we 
marginalized over A (or R) from 10 to 100 Mpc/h. 

RESULTS AND CONCLUSIONS 

Figure shows the Yukawa-like model's parameter 
constraints for the 2dFGRS and SDSS data, and figure 
(0 shows the PL model's constraints. The "ringing" in 
the d function in the PL model, creates "wiggles" in the 
large scale structure power spectrum which are responsi- 
ble for the geometry of the iso-likelihood contours in fig 3. 
The left plots show the likelihood contours for one-sigma 
marginalized, one-sigma joint, and two-sigma joint con- 
fidence levels for one survey in A, a (or R,e) space. The 
center plots show the difference in chi square (A% 2 ) as a 
function of a or e, marginalized over A or R. We used 
data up to k ~ 0.15 h/Mpc for both surveys. The results 
are consistent with each other and with a strict inverse- 
square law. Since the constraints are weak and the Tay- 
lor expansion analysis breaks down at large values of \a\ 
and |e|, we extrapolated the error bars from the graphs. 
We find a = 0.025±^, e = 0M±H for 2dFGRS and 
a = -OMtoi, e = 0A5±{\i for SDSS, at the one-sigma 
level. Non-linearities become increasingly important as k 
increases and might introduce systematic errors, but the 
signal-to-noise for the power spectrum increases with in- 
creasing k. Therefore we explore the dependence of the 
result as a function of the maximum k included in the 
analysis (/c max ) on the right panel of the figures. 

Since our analysis probes deviations in the shape of 
today's power spectrum, the strongest constraints on a 
and e occur at scales where the slope of d (with a = 1) is 
greatest. This typically corresponds to A and R towards 
the smallest scales of the data, so the smallest comoving 
scales probed have just reached the physical scale of the 
force-law deviation today, and the larger comoving scales 
have passed it. (A deviation near a physical scale A or R 



has only affected comoving scales that have grown larger 
than A or R.) The weakness of our constraints arises pri- 
marily from the need to marginalize over a wide range of 
power spectrum amplitudes and the spectral slope. Fu- 
ture weak lensing surveys will measure the power spec- 
trum for dark matter directly. If a(t) can also be mea- 
sured independently, this data will eliminate the need to 
marginalize over the amplitude and shrink our error bars 
for the PL model. We forecast that this would roughly 
halve the error on e. The constraint does not improve 
significantly for the Yukawa-like model. 
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Justification of this assumption, used in all analyses of 
large-scale structure, relies on the statistical properties 
of initially gaussian random fields [2rll27|. 
We neglect the decaying solution. 
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where H and A are functions of a and primes denote 



